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EVALUATION  OF  AN  APPLICATION  OF  ADJOINT  METHODS 
TO  YELLOW  SEA  MODELING 


1.  INTRODUCTION 

Adjoint  methods,  used  for  over  a  decade  in  the  field  of  meteorology,  have  provided  a  powerful 
tool  to  ocean  modelers  as  field  observations  are  more  readily  available  and  data  assimilation  be¬ 
comes  more  computationally  possible.  The  evolution  and  review  of  some  applications  of  adjoints 
and  other  assimilative  methods  can  be  found  in  Ghil  and  Manalotte-Rizzoli  (1991)  and  Ghil  et 
al.  (1997). 

Several  recent  studies  have  applied  adjoint  methods  to  infer  open  boundary  forcings  from 
assimilated  data  (Seiler  1993;  Lardner  1993;  Thompson  and  Griffin  1998).  The  incremental 
approach  (Thompson  et  al.  1998;  Courtier  et  al.  1994;  Thompson  et  al.  2000)  involves  a 
straightforward  modification  of  the  Gauss-Newton  algorithm  commonly  used  for  iteratively  solv¬ 
ing  nonlinear  regression  problems.  The  adjoint  of  a  dynamically  simple  forward  model  is  calcu¬ 
lated  and  the  cost  function  associated  with  point  data  is  minimized  given  station  elevation  and 
velocity  information.  Further  iterated  with  a  dynamically  complex  ocean  model,  this  approach 
can  provide  even  more  accurate  ocean  predictions  relative  to  the  observed  state. 

Adjoint  methods  that  derive  open  boundary  forcing  can  also  facilitate  limited  area  modeling. 
Limited  area  domains  have  the  advantages  of  reducing  computation  time  or  allowing  even  greater 
refinement  in  nearshore  waters.  Traditionally,  though,  this  approach  has  its  shortcomings  in  both 
theoretical  and  practical  aspects.  Blain  et  al.  (1994),  in  a  study  of  the  effect  of  domain  size  on 
computed  storm  surge  response,  found  that  domain  size  and  the  placement  of  the  open  ocean 
boundaries  have  a  direct  effect  on  the  resulting  solution.  Dynamics  generated  within  a  larger 
area  domain  often  cannot  be  captured  within  a  limited  area  domain  or  specified  through  the  open 
boundary  forcing  of  limited  area  domains.  A  practical  advantage  of  large  area  domains  has  been 
that  their  boundaries  extend  far  enough  offshore  that  boundary  forcing  can  be  interpolated  from 
larger  global  solutions.  Open  boundaries  of  limited  area  domains  lie  well  nearshore  of  “proven” 
global  solutions;  thus,  finding  appropriate  open  boundary  forcing  can  become  problematic.  Ad¬ 
joint  methods  may  surpass  such  limitations  to  provide  appropriate  boundary  forcing  for  limited 
area  domains. 

Manuscript  approved  September  25,  2000. 
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For  this  study,  the  complex  forward  model  is  defined  by  a  nonlinear  finite  element  model, 
ADCIRC.  Perhaps  the  greatest  advantage  of  finite  element  discretizations  over  more  traditional 
finite  difference  methods  is  the  ability  to  vary  resolution  over  the  domain,  more  accurately  repre¬ 
senting  the  coastal  boundary  and  increasing  resolution  in  shallower  waters.  In  this  manner,  use 
of  complicated  nested  or  multiresolution  finite  difference  grids  can  be  avoided  while  improving 
the  solution  in  nearshore  regions  and  greatly  reducing  overall  computation  time.  An  available 
adjoint  model  system  of  Thompson  (Griffin  and  Thompson  1996;  Thompson  and  Griffin  1998; 
Pistek  et  al.  1998)  is  selected  as  the  data-assimilative  component. 

The  question  arises  as  to  whether  an  advantage  can  be  gained  by  applying  the  incremental 
approach  to  a  complex  forward  model  whose  adjoint  model  system  is  not  strictly  an  adjoint  of 
the  forward  model.  Certainly,  success  of  this  approach  would  be  appealing,  eliminating  the  need 
to  develop  a  separate  adjoint  for  each  forward  model.  In  this  study,  the  possible  applications 
of  finite  difference-based  adjoint  methods  in  conjunction  with  finite  element-based  models  are 
investigated  in  the  Yellow  Sea.  Specifically,  feasibility  of  the  incremental  approach  is  explored, 
using  the  nonlinear  finite  element  model  ADCIRC  as  the  complex  forward  model  and  the  linear 
finite  difference  adjoint  model  of  Thompson  (ADJOINT)  in  the  linear  inverse  model  loop  with 
data  assimilation.  Also  examined  is  the  degree  to  which  linear  adjoint  methods  can  provide 
more  accurate  open  boundary  forcing  for  nonlinear  forward  models  applied  over  limited  area 
domains  external  to  the  data  assimilative  loop.  After  discussing  the  specifics  of  the  model  system, 
the  report  presents  the  finite  element  grids  and  observational  station  data  used  in  the  model 
system.  The  later  sections  describe  and  analyze  the  implementation  of  both  non-assimilative 
and  assimilative  experiments  and  present  conclusions. 

2.  MODEL  SYSTEM 

The  incremental  modeling  approach  involves  two  components.  First,  a  dynamically  simple 
model  and  its  adjoint,  in  this  case,  the  ADJOINT  model  of  Thompson  (Griffin  and  Thompson 
1996;  Thompson  and  Griffin  1998;  Pistek  et  al.  1998),  is  used  to  assimilate  data  and  com¬ 
pute  optimal  boundary  forcings  constrained  by  both  the  dynamical  equations  and  observations. 
Secondly,  a  model  which  includes  more  complex  nonlinear  dynamics,  ADCIRC  (Luettich  et  al. 
1992;  Kolar  et  al.  1994a, b;  Westerink  et  al.  1994a, b),  serves  as  the  complex  forward  model. 
The  forward  model  uses  open  boundary  forcing  derived  from  the  solution  of  linear  ADJOINT 
model.  Prediction  errors  in  the  complex  forward  model  solution  at  observational  stations  axe 
then  assimilated  into  the  linear  model  system  to  find  the  increment  by  which  the  controls,  in  this 
case,  forcing  on  the  complex  model  open  boundary,  are  adjusted.  This  process  involves  repeated 
runs  of  the  linear  model  and  its  adjoint  to  iteratively  minimize  the  modified  cost  function,  effec¬ 
tively  giving  a  pair  of  nested  iterations  (Thompson  et  al.  2000).  The  modified  cost  function  is 
essentially  the  difference  between  the  complex  model  prediction  and  the  data  at  the  observational 
points.  This  approach  is  computationally  efficient  and  simple,  but  has  limitations.  In  particular, 
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if  the  simple  model  dynamics  are  not  sufficiently  close  to  the  Jacobian  of  the  complex  model, 
there  is  no  guarantee  of  convergence  to  the  true  state,  or  that  the  solution  will  converge  at  all 
(Thompson  et  al.  1998). 

2.1  ADCIRC 

The  complex  forward  model  computations  are  performed  using  the  finite  element-based 
ADCIRC-2DDI,  the  depth-integrated  version  of  a  set  of  two-  and  three-dimensional  fully  nonlin¬ 
ear  hydrodynamic  codes.  ADCIRC-2DDI  solves  the  two-dimensional,  depth-integrated  shallow 
water  equations,  subject  to  the  hydrostatic  pressure  and  Boussinesq  approximations.  For  the 
applications  discussed  here,  the  baroclinic  terms  in  the  momentum  conservation  equations  are 
neglected,  leading  to  the  following  set  of  governing  equations  in  the  primitive,  nonconservative 
form,  using  spherical  coordinates: 


SC  1  dUH  1  dVH  _  VH  tan<f> 
dt  Rcoscj)  dX  +  R  dcj>  R 
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where  t  is  time,  A,  cj)  are  degrees  longitude  (east  of  Greenwich  positive)  and  degrees  latitude 
(north  of  the  equator  positive),  C  is  free  surface  elevation  (relative  to  the  geoid),  U,V  are  depth 
integrated  horizontal  velocities,  and  H  is  the  total  water  column  thickness,  f 2  is  the  angular 
speed  of  the  Earth  (7.29212xl0~5  rad/s),  R  the  radius  of  the  Earth,  /  (=  2f!sin  cj))  is  the  Coriolis 
parameter  determined  by  latitude,  g  is  gravity,  and  (7?  +  T)  represents  the  Newtonian  tidal 
potential,  earth  tide,  self  attraction,  and  load  tide  term.  Lastly,  p0  is  the  reference  density  of 
water.  Standard  quadratic  parameterizations  for  the  bottom  stresses  are  used: 


TbX  =  Cf(U2  +  V2)1?2U 


(4) 
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r b<p  =  Cf(U2  +  V2)1/2V, 


(5) 


where  Cj  is  a  dimensionless  friction  coefficient. 

The  Generalized  Wave  Continuity  Equation  is  obtained  by  reformulating  Eqs.  (1)  through 
(3)  and  is  discretized  in  space  using  the  finite  element  method  (Lynch  and  Gray  1979;  Lynch 
1983;  Kinnmark  1984;  Westerink  and  Gray  1991)  and,  in  time,  using  finite  difference  techniques. 
Further  details  of  the  ADCIRC  family  of  codes  are  presented  by  Luettich  et  al.  (1992),  Kolar  et 
al.  (1994a, b),  and  Westerink  et  al.  (1994a, b). 

2.2  ADJOINT 

In  the  application  of  adjoint  methods,  the  adjoint  of  the  discrete  linear  simple  model  equa¬ 
tions  is  used  to  compute  the  directions  in  which  unknown  model  parameters  are  adjusted  to 
achieve  a  better  fit  of  the  modeled  dynamics  to  observed  values.  The  equations  of  motion  used 
by  the  linear  ADJOINT  forward  model  are  expressed  in  Cartesian  coordinates  as 


lti~fV  +  9dx  ' 

^  +  fU  +  9fy  =  -T{B'V) 

ac  dim  ■  d{VH) 

dt  dx  dy 


(6) 

(7) 

(8) 


where  the  bottom  stresses  and  are  linear  in  U  and  V)  respectively,  with  frictional 

coefficient  k.  Due  to  widely  varying  bottom  topography  (mud  flats  to  rocky  terrain)  in  the  Yellow 
Sea,  the  linear  bottom  friction  coefficient  k  is  spatially  varied  according  to  five  separate  regimes 
over  the  domain  (Fig.  1).  The  open  boundary  is  forced  with  the  following  condition: 


Un  — 


(9) 


where  h  is  the  bathymetric  water  depth  and  the  time-varying  boundary  forcing  £  is  considered 
unknown  in  the  linear  forward  model  solution.  This  type  of  boundary  condition  allows  transient 
waves  generated  inside  the  model  domain  to  radiate  away  with  minimal  reflection  (Thompson  et 
al.  2000). 


The  unknown  tidal  amplitudes  and  phases  at  the  open  boundary  (  are  estimated  from  as¬ 
similated  interior  station  data  (Fig.  2),  and  the  solution  of  the  shallow  water  equations  is  posed 
as  a  boundary  value  problem  in  the  frequency  domain.  The  estimation  procedure  is  based  on  a 
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Fig,  i  —  Dimensionless  linear  drag  coefficient  tt  used  by  the  linear  ADJOINT  model  in  the  Yellow  Sea 
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Fig.  2  —  A  map  of  elevation  stations  and  ADCP  locations  used  for  assimilation  and  comparison  to  model  results 

generalized  least-squares  regression  (Griffin  and  Thompson  1996;  Pistek  et  al.  1998;  Thompson 
and  Griffin  1998). 

3.  APPLICATION  TO  THE  YELLOW  SEA 

The  incremental  approach  is  applied  here  to  tidal  dynamics  in  the  Yellow  Sea,  the  semi- 
enclosed  sea  between  China’s  east  coast  and  Korea  (Fig.  3).  This  marginal  sea  is  one  of  the 
largest  shallow  areas  of  continental  shelf  in  the  world.  In  addition,  there  exist  sufficient  observed 
data  to  effectively  use  adjoint  methods  in  the  region. 
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The  coastline  of  the  region  is  quite  complex,  which  requires  increasing  resolution  nearshore 
for  proper  coastal  representation.  For  this  reason,  applications  of  a  complex  model  over  limited 
domains  using  open  boundary  information  obtained  from  the  incremental  approach  applied  over 
a  larger  grid  are  explored.  Given  proper  information  at  the  open  boundary  of  the  limited  area 
domain,  the  use  of  limited  area  domain  models  can  allow  even  greater  refinement  nearshore 
without  sacrificing  computational  time. 

3.1  Station  Data 

Tidal  elevations  are  obtained  from  the  International  Hydrographic  Office  (IHO)  tidal  sta¬ 
tions,  select  points  along  TOPEX-POSEIDON  altimeter  tracks,  station  data  from  L.  Kantha 
(University  of  Colorado),  and  three  ADCP  records  collected  by  the  U.S.  Naval  Oceanographic 
Office.  Figure  2  is  a  map  showing  the  geographical  locations  and  source  of  data.  Data  at  IHO 
stations  located  beyond  the  boundaries  of  the  finite  difference  grid  used  by  the  linear  ADJOINT 
model,  for  example,  up  rivers  or  inside  harbors,  are  discarded  from  the  dataset  and  are  not  shown 
in  Fig.  2. 

All  comparisons  are  made  for  the  month  of  September  1995,  an  interval  when  winds  were 
weak  and  the  observations  are  dominated  by  tides  (Perkins  and  Pistek  1998).  Time  series  from 
all  stations  were  vertically  integrated  to  make  them  compatible  with  the  assimilation  dynamics, 
then  detided  and  reconstructed  to  a  month-long,  half-hour  dataset  comprising  the  four  main  tidal 
constituents  (Oi,  A; .  M2  ■  and  S2 ) •  This  later  time  series  forms  the  basis  for  comparison  with 
model  solutions. 

3.2  Grids 

Current  NAVOCEANO  tidal  predictions  in  the  Yellow  Sea  are  obtained  from  a  tidal  database 
created  by  the  Army  Corps  of  Engineers.  The  database  is  derived  from  an  ADCIRC  model 
solution  computed  on  a  finite  element  grid  that  is  known  to  under-resolve  the  speed  of  a  gravity 
wave  using  a  dimensionless  wavelength  criterion  (LeProvost  and  Vincent  1986;  Westerink  et  al. 
1994b, c).  This  fareast  grid  (Fig.  4),  with  its  open  boundary  well  off  the  continental  shelf  into  the 
Pacific  Ocean,  was  a  first  attempt  at  creating  a  finite  element  grid  for  applications  in  the  Yellow 
Sea  and  Sea  of  Japan  with  finite  element  hydrodynamic  models  such  as  ADCIRC,  FUNDY,  and 
QUODDY. 

The  ADCIRC  solution  over  the  fareast  grid  shows  mediocre  agreement  with  station  data 
(Fig.  5).  The  solution  at  TOPEX  stations  in  the  interior  of  the  domain  exhibits  poor  agreement 
with  data  (i.e. ,  errors  between  10  and  30  cm).  Areas  where  the  largest  disagreement  with  data 
can  be  found  are  located  off  the  Chinese  coast  near  the  entrance  to  the  Yellow  Sea  and  at  tidal 
stations  near  the  Kyunggi  and  Seahan  Bays  off  the  Korean  coast,  where  errors  exceed  60  cm  in 
root  mean  square  (RMS)  error.  This  deterioration  of  the  solution  along  the  path  of  the  Kelvin 
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Fig.  5  —  Time  series  tidal  elevation  error  by  station  for  the  fareast  grid.  RMS  error  in  meters. 


wave  (Jacobs  et  al.  1998)  from  the  southeast  into  the  interior  and  its  exit  from  the  basin  in  the 
southwest  indicates  a  problem  with  the  grid  itself  rather  than  with  the  solution  method. 


As  a  packet  of  waves  travels  into  progressively  shallower  water,  its  group  velocity  decreases, 
and  greater  resolution  is  necessary  to  resolve  the  waves.  The  fareast  grid  is  inherently  destructive 
in  that  it  does  not  properly  resolve  the  M2  tidal  signal  as  it  propagates  across  the  shelf  and  into 
nearshore  regions  within  the  domain  (Fig.  5).  To  correct  this  deficiency,  a  new  family  of  grids 
for  the  Yellow  Sea-Sea  of  Japan  is  created  to  replace  the  fareast  grid. 

The  yessoj  grid  (Fig.  4),  like  the  fareast  grid,  extends  radially  outward  from  the  Yellow  Sea, 
encompassing  the  East  China,  Philippine,  and  Japan  Seas  and  extending  into  the  deep  water 
of  the  Pacific  Ocean.  The  open  boundary  stretches  from  just  north  of  Hong  Kong,  through  the 
Philippine  Sea  and  deep  waters  of  the  Pacific  Ocean,  to  just  north  of  Tokyo,  Japan  (see  Fig.  3). 
Bathymetry  ranges  from  an  imposed  minimum  of  2  m  in  the  shallowest  waters  of  the  Bohai  Sea  to 
over  8000  m  in  the  Western  Pacific  Ocean  at  the  easternmost  section  of  the  grid.  The  resolution 
in  nearshore  regions,  especially  in  the  Bohai  Sea,  is  fine  tuned  (Fig.  6)  to  resolve  a  range  of  tidal 
frequencies,  increasing  nearshore  resolution  to  unmask  smaller  scale  features  within  the  flow. 
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Resolution  increases  dramatically  at  the  shelf  break,  around  islands  in  the  southwestern  part  of 
the  domain,  and  around  Japan.  On  the  shelf  of  the  Yellow  Sea,  resolution  of  the  yessoj  grid 
(Fig.  6a)  increases  in  the  Bohai  Sea  and  very  close  to  land,  but  decreases  slightly  in  the  center 
of  the  Yellow  Sea,  where  gravity  wave  speed  permits  slightly  lower  resolution  than  that  specified 
by  the  fareast  grid  (Fig.  6b). 

The  bathymetric  depths  for  the  Bohai,  Yellow,  and  East  China  Seas  (depths  less  than  200 
m)  are  determined  from  a  high  resolution  data  set  distributed  by  C.  Horton  at  NAVOCEANO 
in  1995  (Fig.  3).  These  data  originated  from  South  Korean  sources  and  were  compiled  and 
processed  by  L.  H.  Kantha  at  the  University  of  Colorado.  The  spatial  resolution  for  this  data 
is  5  minutes  in  longitude  and  4  minutes  in  latitude.  For  locations  where  depths  exceed  200 
m,  bathymetry  is  derived  from  the  NAVOCEANO  Digital  Bathymetric  Data  Base  -  5  minute 
resolution  (DBDB-5).  Additionally,  an  updated  bathymetry  was  included  in  the  Yellow  Sea  north 
of  the  34°  N  line  of  latitude.  Adjustments  to  the  bathymetry  are  made  based  on  preliminary 
ADJOINT  experiments.  These  modifications  are  included  for  compatibility  between  the  forward 
and  adjoint  model  bathymetry. 

From  the  yessoj  grid,  a  family  of  grids  (Fig.  4)  is  created  with  open  boundaries  positioned 
successively  farther  into  the  Yellow  Sea  on  the  continental  shelf.  All  domains  used  here  extend 
outward  from  the  Bohai  Sea  and  cover  areas  of  decreasing  size.  By  retaining  the  same  resolution 
throughout  the  domain  and  changing  only  the  placement  of  the  open  boundary,  the  effect  of  large 
vs  limited  area  domains  can  be  investigated  in  conjunction  with  ADJOINT  model  experiments. 

4.  MODEL  EXPERIMENTS 

4.1  Non-assimilative  Experiments 

The  initial  yessoj  solution  is  forced  by  tidal  potential  forcing  and  specified  elevations  on  the 
open  boundary  derived  from  eight  tidal  frequencies  taken  from  the  Grenoble  global  tidal  database 
(Le  Provost  et  al.  1994).  Though  the  results  were  analyzed  with  respect  to  only  four  primary 
constituents,  the  M2,  S2,  0\  and  K\,  the  K2,  N2,  Pi,  and  Qi  tidal  frequencies  are  included  in 
open  boundary  forcing  for  completeness. 

Figure  7  shows  the  cotidal  chart  for  M2  dynamics.  Two  amphidromes  are  located  off  the 
Chinese  coast,  and  another  just  offshore  from  the  entrance  into  the  Bohai  Sea.  M2  elevation 
increases  as  the  tidal  wave  propagates  onto  the  shelf,  up  the  Korean  coast,  and  into  Kyunggi 
Bay  near  Seoul  (Fig.  3).  Amplification  of  the  M2  tide  is  also  pronounced  in  the  Tsushima  Strait 
separating  Taiwan  and  mainland  China.  Current  residuals  through  the  domain  are  generally 
small  (3  to  6  cm/s),  except  near  islands  in  the  East  China  and  Japan  Seas,  where  they  range 
from  10  to  25  cm/s  in  magnitude.  These  larger  residual  currents  may  be  an  indication  of  under¬ 
resolution  of  the  mesh  in  these  regions  of  steep  topographic  gradients. 
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Fig.  8  —  Difference  (y es s oj- f areas t)  in  mean  absolute  error  by  station.  Magnitude  of  the  difference  is  color-coded  in 

meters. 

Comparison  with  the  observational  station  data  shows  significant  improvement  of  mean  ab¬ 
solute  error  over  the  fareast  solution  at  all  80  stations  analyzed  (Fig.  8),  ranging  from  a  5  cm 
improvement  through  the  interior  of  the  domain  to  as  much  as  50  cm  at  the  tidal  stations  near 
the  coast.  The  greatest  improvement  was  seen  at  the  tidal  stations  nearest  the  coastlines  in 
shallow  water,  where  the  solution  improved  on  the  order  of  50%  to  80%,  from  70  cm  to  20  cm  in 
mean  absolute  error. 


The  finite  element  grid  adjoint  (Fig.  4)  is  created  to  match  exactly  the  boundaries  of  the 
finite  difference  grid  used  with  the  linear  ADJOINT  model.  The  initial  open  boundary  of  this 
grid  stretches  from  the  coast  of  China  at  Yancheng,  well  inside  the  shelf-break,  at  approximately 
33.25°  N  latitude.  Interpolation  from  the  global  tidal  database  used  previously  (Grenoble)  is 
clearly  inappropriate  for  a  boundary  located  this  far  onto  the  continental  shelf.  Thus,  the  yessoj 
solution  is  interpolated  to  the  adjoint  open  boundary  nodes  and  used  as  elevation  forcing  on  the 
boundary. 
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Fig.  9  —  The  open  boundaries  of  the  initial  (blue)  and  final  (black)  versions  of  the  adjoint  grid 

The  combination  of  the  boundary  placement  through  Cheju  Island,  and  a  90°  corner  in  the 
open  boundary  of  the  new  grid  proved  problematic  (blue  line  in  Fig.  9).  Circulation  around 
Cheju  Island  which  had  developed  in  the  yessoj  solution  now  used  as  forcing  is  constrained  by 
the  essential  elevation  boundary  condition  as  applied  in  the  adjoint  model.  A  nonphysical  com¬ 
pensating  flow  in  the  tidal  residual  currents  emerges,  and  exhibits  signs  of  numerical  instability 
because  of  the  sharp  corner  and  ill-posed  boundary  specifications. 

A  series  of  experiments  having  varied  specification  of  a  radiation  boundary  condition  along 
the  open  ocean  boundaries  of  the  adjoint  grid  aim  to  recapture  a  dynamically  consistent  residual 
flow  near  Cheju  Island.  Instead  of  specifying  elevation  on  the  open  boundary,  the  boundary  is 
treated  as  a  natural  boundary  allowing  wave  radiation  normal  to  the  boundary  (Luettich  et  al. 
1992;  Westerink  et  al.  1994a).  This  radiation  condition  is  applied  by  specifying  the  contribution 
to  the  normal  boundary  flux  integral  in  the  continuity  equation,  a  formulation  that  is  unique  to 
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the  finite  element  Generalized  Wave  Continuity  Equation-based  models  (Luettich  et  al.  1992; 
Kolar  et  al.  1994a, b;  Westerink  et  al.  1994b).  In  this  formulation,  there  is  no  constraint  on 
velocity  (normal  or  tangential)  in  the  momentum  equations;  thus,  waves  are  allowed  to  propagate 
freely  out  of  the  domain.  The  flow  normal  to  the  boundary  is  computed  using  a  Sommerfeld 
radiation  condition  (Roed  and  Smedstad  1984).  This  boundary  condition  correctly  satisfies  the 
flux  balance  in  a  global  sense,  but  only  satisfies  the  normal  flow  at  each  boundary  node  in  the 
limit  of  infinite  resolution  (Westerink  et  al.  1994b). 

The  flux  specified  as  part  of  the  radiation  boundary  condition  implementation  is  interpolated 
from  the  yessoj  solution  in  the  same  manner  as  for  specified  elevations.  Four  radiation  boundary 
regimes  are  considered:  a)  contiguous  implementation  on  the  open  boundaries  north  and  east 
of  Cheju  island,  b)  implementation  separately  along  the  north  and  east  segments  of  these  open 
boundaries,  c)  implementation  of  the  radiation  boundary  condition  on  the  north  segment  and 
specifying  elevation  on  the  east  segment  (and  vice  versa),  and  d)  implementation  of  a  contiguous 
radiation  boundary  condition  as  in  a),  but  now  “tacking  down”  the  solution  at  the  juncture  of  the 
north  and  east  boundary  segments  with  five  specified  elevation  boundary  nodes.  Despite  these 
various  configurations  of  the  radiation  boundary  conditions,  the  90°  corner  in  the  grid  remained 
problematic  in  its  generation  of  a  numerically  spurious  solution.  As  a  final  solution,  the  grid  was 
modified,  moving  the  eastern  boundary  westward  such  that  the  boundary  now  extends  from  the 
Korean  coast  directly  south  to  Cheju  Island  (black  line  in  Fig.  9). 

As  mentioned  previously,  bottom  topography  in  the  Yellow  Sea  ranges  from  mud  flats  in 
the  Bohai  Sea  to  rocky,  jagged  terrain  off  the  South  Korean  coast.  Recall  that  the  linear  AD¬ 
JOINT  model  uses  spatially  varying  coefficients  for  the  linear  friction  formulation  over  the  domain 
(Fig.  1).  Dimensionless  values  range  from  0.01  in  the  Liaodong  Bay  in  the  northeastern  Bohai 
Sea  to  0.18  off  the  South  Korean  coast  in  the  east.  To  investigate  the  possibility  of  including 
spatially  variable  friction  coefficients  for  use  with  the  quadratic  bottom  stress  formulation  in 
ADCIRC,  a  series  of  model  runs  are  conducted  with  a  spatially  constant  friction  coefficient,  Cf, 
set  to  0.0015,  0.002,  and  0.0025.  The  results  from  these  three  runs  are  compared  to  station  data. 

RMS  error  is  calculated  by  station,  and  analyzed  with  respect  to  station  type  and  geograph¬ 
ical  location.  In  general,  it  is  found  that  data  at  most  stations  agree  best  with  model  results 
using  frictional  values  of  Cf  —  0.002  and  0.0025.  The  exception  is  in  the  middle  of  the  basin 
(Fig.  10),  where  the  friction  coefficient  of  0.0015  resulted  in  smaller  RMS  errors.  It  is  interesting 
to  note  that  the  observations  at  stations  in  the  Bohai  Sea  agree  best  with  the  friction  values 
of  0.002  and  0.0025.  Minimum  errors  at  TOPEX  track  stations  are  recorded  using  the  lowest 
friction  coefficient,  Cf  =0.0015. 
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Fig.  10  — ■  Minimum  error  by  station  for  varying  Cf.  Color  of  the  error  shape  represents  the  magnitude  of  the  minimum  root  mean 
square  error  in  meters.  Shape  represents  which  frictional  value  produced  that  minimum. 
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Prom  this  set  of  experiments,  it  becomes  clear  that  specifying  complicated  spatially  variable 
drag  coefficients  over  the  majority  of  the  domain  provides  relatively  little  improvement.  The 
nearshore  regions  off  the  Korean  coast  post  the  largest  RMS  errors  yet  appear  relatively  in¬ 
sensitive  to  the  frictional  coefficient.  For  all  remaining  ADCIRC  model  simulations,  a  constant 
frictional  coefficient  of  0.0025  is  used. 

4.2  Assimilative  Experiments 

With  the  open  boundary  of  the  adjoint  grid  coincident  with  the  open  boundary  of  the  finite 
difference  ADJOINT  model,  the  complex  forward  model  ADCIRC  is  forced  using  open  boundary 
values  generated  by  the  solution  of  the  simple  linear  ADJOINT  model.  Model-data  misfits 
between  the  forward  model  predictions  and  observed  elevation  and  velocity  time  series  at  data 
stations  are  assimilated  by  the  ADJOINT  model.  Optimal  boundary  conditions  are  then  derived 
by  minimizing  this  error  between  the  forward  model  and  the  data. 

Flux  from  the  ADJOINT  model  solution  is  interpolated  to  the  open  boundary  of  the  adjoint 
grid  for  each  of  the  four  constituents  previously  analyzed,  the  M2,  52,  0i,  and  K\.  The  difference 
in  resolution  along  the  boundary  of  the  two  grids  (40  km  vs  4-20  km)  results  in  discrepancies 
in  the  land/sea  interface.  As  such,  it  is  necessary  to  extrapolate  the  fluxes  at  the  eastern  and 
westernmost  open  boundary  nodes  of  the  ADJOINT  model  to  the  adjoint  grid  nodes  that  extend 
to  the  mainland  using  a  band-limited  reconstruction  technique.  In  this  approach,  the  magnitude 
of  the  flux  is  tapered  sinusoidally  towards  a  mean  flux  (computed  from  the  linear  ADJOINT 
model  solution)  at  the  mainland.  Despite  the  large  differences  in  resolution,  flux  from  the  linear 
ADJOINT  model  solution  is  sufficiently  smooth  for  interpolation  to  interior  grid  points  along  the 
adjoint  finite  element  grid  open  boundary. 

Cast  as  the  complex  forward  model  in  iteration  with  the  simple  ADJOINT  model,  ADCIRC 
is  then  run,  forced  by  this  specified  normal  flux  on  the  open  boundary.  Within  ADCIRC,  the 
normal  flux  is  represented  by  essential  boundary  condition  with  no  constraint  on  tangential 
flow  at  the  open  boundary  nodes.  This  boundary  condition  implemention  requires  specifying 
a  nonzero  contribution  to  the  normal  boundary  flux  integral  in  the  continuity  equation  and  a 
nonzero  normal  velocity  in  the  momentum  equations.  A  flux  balance  is  achieved  in  a  global 
sense  between  the  two  models  while  satisfying  the  need  for  a  normal  flux  at  each  boundary 
node.  Upon  completion  of  an  ADCIRC  model  simulation  on  the  adjoint  grid,  a  time  series  at 
each  station  is  created  from  the  ADCIRC  elevation  and  velocity  solutions.  The  linear  ADJOINT 
model  is  then  run  with  these  time  series  representing  deviations  from  the  observational  data  used 
previously.  This  latter  ADJOINT  model  solution  at  the  open  boundary  is  then  reinterpolated 
and  extrapolated  to  the  complex  model  adjoint  grid  boundary  and  the  iteration  begins  again. 

In  this  application  to  the  Yellow  Sea,  each  successive  iteration  yields  model-data  misfits 
between  the  forward  model  and  the  observations  that  grow  in  magnitude,  indicating  divergence 
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of  the  solutions.  To  test  the  role  of  the  band-limited  extrapolation  scheme  in  this  divergence, 
boundary  values  at  grid  nodes  in  the  extrapolation  region  are  assigned  velocities  interpolated 
from  the  yessoj  solution  (Fig.  11).  This  scheme  improves  the  forward  model  solution  slightly 
with  respect  to  RMS  errors,  but  when  further  iterated  with  the  ADJOINT  model,  the  solutions 
continue  to  diverge. 

Given  the  experiences  thus  far  and  knowing  that  the  incremental  data  assimilative  approach 
is  limited  by  the  proximity  of  the  simple  model  dynamics  to  the  Jacobian  of  the  complex  model, 
a  simple  twin  experiment  is  designed  to  determine  whether  convergence  of  this  forward- inverse 
model  system  is  possible.  An  ADCIRC  solution  is  obtained  for  a  simple  M2  flux  forcing  across  the 
open  boundary.  Assuming  that  the  ADCIRC  model  predictions  exactly  match  the  observations 
at  the  80  locations  previously  identified,  the  ADJOINT  model  should  recover  the  M2  flux  forcing 
along  its  open  boundary.  However,  the  linear  ADJOINT  model  is  unable  to  recover  this  simple 
forcing  at  its  open  boundary  with  the  ADCIRC  computed  elevation  and  velocity  time  series 
assimilated  at  80  observation  points  (personal  communication,  J .  Book) . 

This  result  suggests  that  convergence  between  the  two  models  in  this  application  is  unlikely 
and  may  be  the  result  of  known  incompatibilities  between  the  models.  It  is  well  known  that 
a  properly  formulated  ADJOINT  model  is  derived  from  the  linearized  discrete  form  of  the  for¬ 
ward  model  (Bennett,  1992;  Thompson  et  ai,  2000).  This  ensures  that  the  minimization  search 
proceeds  in  a  direction  appropriate  for  determination  of  the  true  minimum.  ADCIRC- 2DDI  is 
based  on  a  Generalized  Wave  Continuity  Equation-based  formulation  of  the  shallow  water  equa¬ 
tions  and  discretized  using  finite  element  methods  (Luettich  et  al.,  1992;  Kolar  et  al.,  1994a, b; 
Westerink  et  al.,  1994a).  In  contrast,  the  ADJOINT  model  is  a  finite  difference  discretization  of 
the  linearized  primitive  shallow  water  equations.  While  the  forward  and  inverse  adjoint  models 
of  the  linear  system  are  well-posed,  the  inverse  of  the  ADJOINT  model  is  not  a  proper  adjoint 
formulation  for  the  complex  forward  model.  Though  this  incompatibility  was  known  at  the  on¬ 
set  of  this  application  to  the  Yellow  Sea,  the  question  remained  whether  this  incompatibility 
precludes  improvements  to  the  forward  model  solution  from  such  an  implementation  of  the  in¬ 
cremental  approach.  An  answer  to  this  question  is  still  unknown  as  there  are  several  other  model 
incompatibilities  that  may  be  affecting  the  success  of  the  methodology. 

Several  differences  between  the  models  complicate  the  issue  of  convergence  using  the  in¬ 
cremental  approach.  Namely,  the  resolutions  and  computational  grids  of  the  complex  forward 
model  and  the  simple  ADJOINT  model  are  disparate.  Computations  for  the  ADCIRC  model 
are  obtained  over  a  highly  resolved  triangular  mesh  (resolution  ranges  from  4  to  20  km  along  the 
open  ocean  boundary).  The  ADJOINT  model  solution  is  relative  to  an  Arakawa  C-grid  having 
40  km  resolution,  almost  l/10th  the  resolution  of  the  complex  forward  model  grid. 

It  should  be  noted  that  for  the  strong  constraint  inverse  problem,  the  linear  AD  J  OINT  model 
exactly  matches  linear  dynamics  in  the  interior  of  the  domain.  At  the  open  boundary,  an  optimal 
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Fig.  11  —  Mi  velocity  corresponding  the  specified  flux  applied  at  the  southern  (a)  and  eastern  (b)  open 
boundaries  of  the  adjoint  grid.  The  extrapolation/interpolation  scheme  is  shown  in  red,  and  yessoj  solution 
substitution/interpolation  in  blue. 
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solution  is  derived  subject  to  the  dynamical  equations  and  a  minimization  of  the  forward  model- 
data  misfit.  Comparing  the  original  yessoj  solution  at  the  adjoint  grid  boundary  locations  to  the 
fluxes  generated  by  the  solution  of  the  linear  ADJOINT  model  system  illustrates  the  difference 
between  the  optimal  and  purely  dynamical  model  boundary  forcings  (Fig.  12).  “Iteration  1” 
velocities  in  Fig.  12  represent  boundary  forcings  derived  from  the  linear  ADJOINT  model  given 
model-data  misfits  obtained  using  ADCIRC  elevations  and  velocities  (as  computed  using  the 
open  boundary  forcing  shown  in  red).  Large  gradients  in  normal  flux  forcing  are  evident,  even 
as  much  as  a  180  degree  phase  change  over  less  than  100  km  in  the  westernmost  section  of  the 
open  boundary,  and  the  magnitude  of  this  gradient  grows  substantially  between  iterations  of  the 
linear  ADJOINT  model  system.  It  is  likely  that  the  coarse  resolution  of  the  ADJOINT  model 
is  a  primary  cause  and  the  conflicting  land/sea  boundary  locations  are  a  secondary  cause  of  the 
boundary  value  discrepancies  depicted  in  Fig.  12.  The  effect  is  divergence  of  the  incremental 
approach  applied  to  the  ADCIRC  and  ADJOINT  models  as  configured. 

Lastly,  the  differences  in  the  implementation  of  the  boundary  conditions  may  also  contribute 
to  the  divergence  of  the  incremental  approach  as  applied  here.  When  forcing  with  specified 
elevation  or  flux  as  an  essential  boundary  condition,  ADCIRC  regards  such  information  at  the 
open  boundary  nodes  as  known.  Thus,  in  the  Generalized  Wave  Continuity  Equation  solution 
matrix,  the  equations  associated  with  those  boundary  nodes  are  eliminated.  In  contrast,  the  open 
boundary  condition  of  the  linear  ADJOINT  model  radiates  about  both  elevation  and  velocity 
at  the  open  boundary  as  calculated  by  the  inverse  model,  using  the  boundary  condition  (9)  to 
minimize  reflection.  This  mismatch  in  boundary  condition  formulation  between  the  complex 
forward  model  and  the  linear  ADJOINT  model  likely  contributes  to  increasing  the  distance 
between  the  Jacobian  of  the  complex  model  and  that  of  the  simple  ADJOINT  model.  According 
to  Thompson  et  al.  (1998),  this  effect  can  be  destructive  to  the  convergence  of  the  solution  using 
the  incremental  approach. 

In  summary,  several  factors  limit  the  success  of  the  data  assimilative  incremental  approach 
as  applied  here.  These  include  incompatibility  of  the  complex  forward  model  and  the  simple 
ADJOINT  model,  disparity  in  resolution  of  the  grids  used  by  the  complex  and  simple  forward 
models,  and  different  boundary  condition  implementations.  The  extent  to  which  each  of  these 
factors  influences  convergence  of  the  incremental  approach  is  yet  to  be  determined.  Certainly, 
this  work  demonstrates  the  need  to  consider  these  issues  of  compatibility  when  implementing  a 
data  assimilative  approach. 

4.3  Limited  Area  Domains 

One  alternative  to  the  incremental  approach  that  includes  a  complex  forward  model  explicitly 
within  the  assimilation  loop  is  to  simply  interpolate  the  optimal  solution  obtained  from  the  linear 
forward/inverse  ADJOINT  model  system  to  the  open  boundary  of  the  complex  forward  model 
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Fig.  12  —  M2  velocities  on  the  southern  (a)  and  eastern  (b)  open  boundaries  are  shown  for  the  adjoint  solution 
computed  with  yessoj  forcings  (blue),  forcings  interpolated  from  the  linear  ADJOINT  model  solution  (red),  and  the 
ADCIRC  solution  “iteration  1”  (cyan). 
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applied  over  a  limited  domain.  In  this  instance,  limited  domain  implies  that  the  nonlinear  forward 
model  open  ocean  boundary  is  contained  within  the  linear  ADJOINT  model  domain. 

In  order  to  test  this  approach  for  obtaining  boundary  forcings  for  the  nonlinear  model,  a 
family  of  limited  area  domains  is  created  from  the  original  yessoj  grid,  with  the  open  boundaries 
positioned  successively  further  into  the  Yellow  Sea  on  the  continental  shelf  (Fig.  13).  Tidal 
elevations  are  interpolated  from  the  linear  ADJOINT  solution  to  the  open  boundaries  of  the 
limited  area  grids  for  each  of  the  four  constituents  previously  analyzed.  The  complex  model, 
ADCIRC,  is  then  executed  over  each  of  these  grids  using  this  interpolated  forcing.  The  results 
are  then  compared  to  the  linear  ADJOINT  model  solutions  that  include  data  assimilation. 

yesbit 

Figure  14  shows  the  difference  in  mean  absolute  error  between  the  ADCIRC  solution  of 
the  yesbit  grid,  produced  with  forcing  interpolated  from  the  original  yessoj  solution,  and  the 
ADCIRC  solution  forced  by  the  interpolated  linear  ADJOINT  model  solution  (denoted  “yesbit 
il”).  Positive  values  indicate  that  greater  errors  occur  in  the  nonassimilative  yesbit  solution. 
Differences  are  small  in  the  interior  of  the  basin  and,  surprisingly,  in  the  Bohai  Sea,  but  increase 
in  magnitude  along  the  eastern  and  western  shores  of  the  Yellow  Sea,  with  greatest  differences 
at  tidal  stations  just  off  the  coast. 

subadj 

The  open  boundary  of  the  subcirdj  grid  is  selected  as  a  straight  line  between  an  altimeter 
station  and  a  Kantha  elevation  station.  Solutions  over  this  finite  element  grid  are  forced  in  a 
similar  manner  as  for  the  yesbit  experiment.  Elevations  at  open  boundary  nodes  are  derived 
by  interpolation  from  the  interior  solution  of  the  linear  ADJOINT  model.  Figure  15  shows  the 
difference  in  error  between  the  solutions.  Again,  positive  differences  represent  greater  error  in 
the  subadj  solution  forced  by  yessoj  elevations.  Differences  in  error  between  the  solutions,  for  the 
most  part,  are  quite  small  (<1.0  cm  RMS).  However,  for  this  grid,  the  original  subadj  solution 
outperforms  the  data  assimilative  (“il”)  solution  that  uses  46  of  the  58  total  stations  that  lie 
in  the  subadj  domain  giving  smaller  RMS  errors  (elevation).  The  exceptions  are  at  two  stations 
that  lie  along  the  open  boundary,  where  errors  of  20  cm  are  recorded. 

subadj2 

A  third  limited  area  domain  is  defined  {subadj 2),  with  its  open  boundary  located  even  further 
into  the  interior  of  the  Yellow  Sea.  Solutions  are  computed  in  an  analogous  manner  as  previous 
limited  area  domains,  and  the  errors  are  analyzed  for  both  simulations.  Fig.  refsub2  shows  a 
comparison  of  RMS  errors  of  the  results. 

As  with  the  comparisons  for  the  subadj  grid  solutions,  use  of  the  assimilative  open  boundary 
forcing  improves  the  solution  significantly  on  the  open  boundary.  However,  unlike  the  subadj 
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Fig.  13 


Limited  area  domain  grids,  plotted  with  the  adjoint  grid  used  in  previous  ADJOINT  experiments 
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grid  solutions,  the  use  of  open  boundary  forcing  from  the  linear  ADJOINT  model  improves  the 
solution  throughout  nearly  all  of  the  domain  (at  the  boundary,  through  the  interior,  and  into 
the  Bohai  Sea).  Near  observation  stations,  the  solution  is  constrained  by  the  assimilation  of 
data  values  (Pistek  et  al.  1998).  The  improvement  in  this  subadj2  solution  at  this  particular 
boundary  may  be  the  result  of  the  proximity  of  the  open  boundary  to  many  locations  where  data 
are  assimilated  into  the  linear  ADJOINT  model. 

The  results  of  these  three  limited  area  domain  experiments  demonstrate  that  the  linear 
ADJOINT  model  solution  is  highly  dependent  upon  the  observational  data  assimilated.  Three 
limited  area  domains  having  three  different  open  ocean  boundary  locations  result  in  complex 
model  solutions  that  are  unpredictable  with  respect  to  their  improvement  or  degradation  in 
comparison  to  measured  data.  The  coarse  resolution  of  the  ADJOINT  model  as  applied  here  is 
a  likely  cause  of  this  data  dependence. 

5.  CONCLUSIONS 

An  adjoint  model  system  comprising  a  dynamically  simple  model  and  its  inverse  in  con¬ 
junction  with  a  more  dynamically  complex  forward  model  has  been  applied  in  the  Yellow  Sea 
to  predict  tidal  elevations  and  currents  throughout  the  region  using  the  incremental  approach. 
A  systematic  series  of  model  experiments  has  shown  that  the  incremental  approach  diverges 
for  the  model  configurations  employed.  Several  factors  contribute  to  this  result,  foremost  being 
incompatibilities  between  ADCIRC  and  the  Thompson  linear  ADJOINT  model,  as  configured. 
The  two  models  solve  different  model  equations,  have  radically  different  discrete  formulations, 
widely  disparate  resolution,  and  differing  open  boundary  forcing  representations.  Limited  do¬ 
main  modeling  using  the  assimilative  model  as  a  source  of  boundary  forcing  is  also  investigated. 
This  approach  indicates  a  strong  dependence  of  the  linear  ADJOINT  model  as  defined  upon  the 
observational  data. 

One  of  the  primary  impedances  to  success  of  the  incremental  approach  as  applied  is  the  widely 
disparate  resolution  of  the  computational  grids  used  by  the  complex  forward  model  and  the  linear 
ADJOINT  model.  The  complex  model  ADCIRC  is  exercised  over  a  mesh  with  resolution  from  4 
to  20  km,  sufficient  to  resolve  the  speed  of  gravity  waves.  In  contrast,  the  linear  (finite  difference) 
ADJOINT  model  has  a  far  coarser  resolution  of  40  km,  which  is  demonstrated  to  be  sufficient  for 
a  linear  model  constrained  by  assimilated  observations  (Perkins  and  Pistek  1998).  This  coarse 
resolution  is  likely  insufficient  without  data  assimilation  and  increases  the  data  dependence  of 
the  linear  model  solution.  Furthermore,  the  differences  in  resolution  create  a  land/sea  boundary 
mismatch  between  forward  and  ADJOINT  models  that  contribute  further  to  divergence  of  the 
method.  Typical  applications  of  adjoint  models  using  the  incremental  approach  maintain  identical 
resolutions  for  the  forward  and  inverse  models. 
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Another  discrepancy  between  the  ADCIRC  and  the  Thompson  model  is  the  theoretical  dif¬ 
ference  in  the  boundary  condition  formulations.  Using  the  incremental  approach  as  applied  here, 
the  linear  ADJOINT  model  is  forced  by  surface  elevations  as  well  as  fluxes,  and  the  open  bound¬ 
ary  condition  is  radiated  about  that  state  imposed  at  the  boundary.  By  contrast,  the  complex 
forward  model,  ADCIRC,  uses  specified  elevation  or  specified  fluxes  along  the  open  boundary. 
No  radiation  of  energy  across  the  open  boundary  takes  place.  In  the  incremental  approach  of 
Thompson  et  al.  (2000),  the  expression  of  the  boundary  forcing  is  an  essential  component  of  the 
optimization  procedure. 

Lastly,  a  known  discrepancy  at  the  start  of  these  experiments  is  the  waiver  on  having  a 
true  adjoint  for  the  complex  forward  model.  The  ADJOINT  model  used  is  not  the  adjoint  of  the 
discrete  linear  equations  of  the  complex  forward  model.  Both  model  equations,  Generalized  Wave 
Continuity  Equation  vs  primitive  Shallow  Water  equations,  and  discretizations,  finite  element  vs 
finite  difference,  are  dissimilar.  These  incompatibilities  are  known  to  be  potentially  problematic. 
This  work,  however,  does  not  prove  conclusively  that  these  factors  are  solely  responsible  for  the 
failure  of  the  incremental  approach.  Rather,  failure  of  the  incremental  approach  as  applied  here 
can  be  ascribed  to  some  combination  of  resolution  disparity,  incompatible  boundary  conditions, 
and  differing  model  formulations. 

In  another  approach,  the  optimal  ADJOINT  model  solution  serves  as  a  source  of  open  bound¬ 
ary  forcing  for  a  limited  area  domain  complex  forward  model  external  to  the  data  assimilation 
loop.  Improvements  in  the  complex  model  performance  as  shown  here  are  unpredictable  due  to 
a  strong  dependence  of  the  ADJOINT  model  solution  on  the  placement  of  the  open  boundary 
and  the  observational  data  assimilated. 

It  is  anticipated  that  an  adjoint  model  system  that  is  derived  from  the  discrete  linear  equa¬ 
tions  associated  with  the  complex  forward  model  would  ensure  a  convergent  solution.  The  forward 
model  specified  by  ADCIRC  and  the  linear  adjoint  model  system  TRUXTON/FUNDY  (Lynch  et 
al.  1998;  Naimie  and  Lynch  2000)  would  be  an  obvious  choice.  Compatible  boundary  conditions 
would  increase  the  likelihood  for  successful  convergence.  Such  a  model  system  could  be  used 
to  further  improve  understanding  of  the  factors  affecting  successful  implementation  of  the  incre¬ 
mental  data  assimilative  approach.  Addition  insight  into  the  tidal  dynamics  and  an  enhanced 
ability  to  forecast  tidal  elevations  in  the  Yellow  Sea  are  also  expected. 
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